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The  resonance  frequency  spectrum  is  first  derived  for  an  m-layered  Goupillaud-type  elastic 
medium  that  is  subjected  to  a  discrete  sinusoidal  forcing  function  that  varies  harmonically 
with  time  at  one  end  and  held  fixed  at  the  other  end.  Analytical  stress  solutions  are  obtained 
from  a  coupled  first-order  system  of  difference  equations  using  z-transform  methods, 
where  the  determinant  of  the  resulting  global  system  matrix  in  the  z-space  |Am|  is  a 
palindromic  polynomial  with  real  coefficients.  The  zeros  of  the  palindromic  polynomial 
are  distinct  and  are  proven  to  lie  on  the  unit  circle  for  1  <  m  <  5  and  for  certain 
classes  of  multilayered  designs  identified  by  tridiagonal  Toeplitz  matrices.  An  important 
result  is  the  physical  interpretation  that  all  the  positive  angles,  coterminal  with  the  angles 
corresponding  to  the  zeros  of  |Am|  on  the  unit  circle,  represent  the  resonance  frequency 
spectrum  for  the  discrete  model.  A  sequence  of  resonance  frequencies  for  the  discrete 
model  appears  to  be  universal  for  all  multilayered  designs  with  an  odd  number  of  layers, 
as  it  is  independent  of  any  design  parameters. 

The  resonance  frequency  results  for  the  discrete  model  are  then  extended  to  describe 
the  resonance  frequency  spectrum  for  the  continuous  model,  where  the  forcing  function 
applied  at  one  end  of  the  strip  is  continuous  and  varies  harmonically  with  time  while 
the  other  end  is  held  fixed.  The  proposed  natural  frequency  spectrum  for  a  free-fixed  m- 
layered  Goupillaud-type  strip  is  confirmed  by  independently  solving  a  simplified  form  of 
the  frequency  equation,  obtained  after  applying  a  transformation  of  the  spatial  variable. 
Our  results  suggest  that  the  natural  frequency  spectrum  depends  on  the  layer  impedance 
ratios  and  it  is  inversely  proportional  to  the  equal  wave  travel  time  for  each  layer. 

The  results  are  used  to  identify  layered  designs  with  a  common  frequency  spectrum  and 
modify  an  existing  design  to  obtain  a  desired  frequency  spectrum.  Other  connections  are 
made  with  previous  stress  optimization  results,  the  Chebyshev  polynomials  of  the  first  and 
second  kind,  as  well  as  the  natural  frequencies  of  a  free-fixed  non-Goupillaud-type  layered 
strip. 

Published  by  Elsevier  B.V. 


1.  Introduction 

The  study  of  natural  vibrations  in  elastic  media  include  the  study  of  resonance,  as  resonance  can  enhance  the  performance 
of  many  sensors  and  devices,  yet  can  devastate  structures  subjected  to  sustained  temporaily-periodic  loading.  Despite  the 
long  history  of  developments  in  the  field,  exact  solutions  for  the  resonance  response  of  multilayered  elastic  media  have  been 
primarily  limited  to  analyses  involving  only  a  few  layers. 
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A  number  of  works  establish  the  natural  modes  of  vibration  in  multilayered  media  [1,2],  and  analogous  problems  for 
n-segmented  strings  [3]  and  n-strings  [4],  Churchill  [5]  uses  Laplace  transform  methods  to  derive  the  transient  resonance 
response  of  a  free-fixed  elastic  bar,  and  notes  the  presence  of  the  product  of  a  temporal  term  and  time-harmonic  function  in 
the  solution  that  becomes  unbounded  at  large  times.  Caviglia  and  Morro  provide  closed-form,  time-domain  expressions  for 
transient  waves  in  a  multilayered  elastic  (possibly  anisotropic)  medium  [6].  They  derive  the  transient  resonance  response 
of  a  single  isotropic  elastic  layer  sandwiched  between  two  halfspaces  that  is  subjected  to  a  temporally  periodic  sawtooth 
function,  and  observe  that  resonance  makes  the  elastic  layer  acoustically  transparent.  A  similar  two-dimensional  problem  is 
studied  by  Kaplunov  and  Krynkin  [7]  who  examine  the  influence  of  layer  stiffness  and  thickness  on  the  resonance  vibrations 
in  a  symmetrically  point-loaded  single  layer  of  elastic  material  enclosed  between  two  half-spaces.  Qiang  et  al.  [8]  study 
the  natural  frequencies  of  anisotropic  multilayers,  and  illustrate  beat  phenomena  in  two  and  three  layer  systems  using 
an  efficient  numerical  eigensolution  scheme  that  is  based  on  semi-analytical  methods.  Exact  expressions  for  the  reflection 
coefficients  for  a  two-dimensional  elastic  layer  overlying  an  elastic  halfspace  are  obtained  by  Fokina  and  Fokin  [9],  but  the 
transient  response  of  the  system  at  resonance  is  not  analyzed.  Graff  presents  a  general  method  [  10]  to  determine  the  natural 
frequencies  of  composite  rods  for  power  ultrasonic  applications  with  a  specific  solution  for  a  two-rod  system.  The  method 
can  be  extended  to  multilayered  systems  but  the  general  problem  of  determining  the  natural  frequencies  with  specific 
geometry  and  material  properties  can  only  be  solved  using  numerical  methods. 

Exact  analytical  expressions  for  the  transient  resonance  response  of  a  multilayered  elastic  medium  are  derivable, 
however,  if  we  specialize  the  medium  to  be  of  Goupillaud-type  [11,12].  The  Goupillaud  specialization  is  often  used  in 
geophysical  applications  to  model  wave  propagation  in  inhomogeneous  media  [13J.  Additional  simplifications  due  to 
transformation  of  the  spatial  variable  to  the  travel-time  variable  have  been  particularly  useful  in  this  work  [14,15].  The 
transformation  has  simplified  the  initial-boundary  value  problem  as  well  as  the  frequency  equation,  allowing  us  to  derive 
analytical  expressions  for  the  natural  frequency  spectrum  in  free-fixed  Goupillaud-type  multilayered  media.  In  addition, 
the  present  treatment  uses  a  global  matrix  method  that  is  attributed  to  Knopoff  [16],  rather  than  the  Thomson-Haskell 
transfer  matrix  formalism  [17,18].  Since  the  recursion  relations  for  the  multilayered  medium  are  written  only  in  terms  of 
the  stresses,  only  half  the  number  of  equations  are  required  relative  to  those  in  classical  global  matrix  methods  [16]. 

Finally,  Bube  and  Burridge  [19]  treat  the  inverse  problem  of  finding  the  coefficients  of  a  first  order  2x2  hyperbolic 
system  related  to  reflection  seismology.  In  order  to  numerically  solve  the  continuous  initial-boundary-value  problem, 
several  difference  schemes  are  applied,  as  discretizations  to  the  corresponding  differential  equations.  The  difference  scheme 
IVp,  given  in  equation  (3.1.11)  pg.  517  of  [19],  appears  to  match  the  recursive  relations  for  the  stress  values  given  in  Section  2. 
The  recursion  relations  and  trigonometric  stress  solutions  are  shown  to  be  exact,  for  discretely  layered  media  of  Goupillaud- 
type.  The  stress  solutions  that  are  derived  herein  using  z-transform  methods,  are  used  to  determine  the  resonance  response 
for  the  discrete  model,  and  subsequently  extended  to  the  continuous  model. 

The  remainder  of  the  paper  is  organized  as  follows.  In  Section  2,  we  derive  stress  solutions  at  resonance  and 
non-resonance  frequencies  for  a  Goupillaud-type  elastic  medium  subjected  to  the  given  boundary  conditions.  This  is 
accomplished  by  generalizing  the  global  system  of  recursion  relations  previously  derived  in  [20]  and  using  z-transform 
methods.  The  determinant  of  the  global  system  matrix  |Am|  is  a  palindromic  polynomial  with  real  coefficients.  We  prove 
that  the  zeros  of  |Am|  are  distinct,  complex  conjugate,  and  lie  on  the  unit  circle  for  1  <  m  <  5  and  for  certain  classes 
of  multilayered  designs  identified  by  tridiagonal  Toeplitz  matrices.  The  zeros  of  the  determinant  allow  us  to  identify  the 
resonance  frequency  spectrum  for  the  discrete  model.  The  resonance  results  derived  for  the  discrete  model  in  Section  2 
are  extended  to  the  continuous  model  in  Section  3,  using  a  linear  transformation  between  the  frequencies  of  the  discrete 
and  continuous  forcing  function.  As  a  result,  we  are  able  to  analytically  describe  the  natural  frequency  spectrum  of  a  free- 
fixed  Goupillaud-type  strip.  An  alternative  derivation  for  the  natural  frequency  spectrum  is  also  provided  in  Section  3 
through  a  generalization  of  Graffs  method  [10].  Common  to  the  analysis  in  Sections  2  and  3  is  the  transformation  and 
simplification  of  the  initial  boundary  value  problem  which  replaces  the  spatial  variable  x  with  the  travel-time  variable  £. 
Various  applications  and  interpretations  are  discussed  in  Section  4.  Section  5  summarizes  the  results,  while  miscellaneous 
references  and  derivations  are  included  in  the  Appendices. 


2.  Discrete  forcing  function  and  resonance  frequencies 

2.3.  Description  of  the  discrete  model 

A  finite  strip  made  of  m-layers  of  homogeneous  elastic  materials  is  considered,  where  the  density  and  elastic  modulus  of 
the  material  occupying  the  ith  layer  are  represented  by  p,  and  £,  respectively,  for  i  =  1 , ,m.  We  relate  to  the  ith  layer,  the 
wave  speed  c,-  =  (Fj/p,)1^2  and  the  characteristic  impedance  p,-  •  q,  for  i  =  1, . . . ,  m.  The  material  properties  and  the  wave 
speed  are  piecewise  constant  functions,  taking  constant  values  in  each  layer.  The  impedance  ratio  between  layers  i  and  i  +  1 
is  represented  by  ctf  =  — ^ —  =  -xSE  ,  for  i  =  1, _ m  —  1,  while  the  transit  time  through  the  strip  is  denoted  by  r. 

ci+l  A'+l  -y/q+lA-l-l 

The  m-layered  strip  is  of  Goupillaud-type,  which  means  equal  wave  travel  time  of  ^  for  each  layer.  The  strip  is  subjected  to 
zero  initial  conditions,  a  discrete  loading  function/(n),  at  one  end,  with  n  >  0  and  step  — ,  and  held  fixed  at  the  other  end, 
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Fig.  1.  Lagrangian  diagram  for  an  elastic  strip  made  of  m-layers  of  equal  wave  travel  time. 


see  Fig.  1.  As  a  result,  our  initial-boundary  value  problem  is 


3  2u(x,  t)  2  3  2u(x,t) 

3 12  = *  3x2  ’ 

a(0,t)  =  E^(0,t)=f(n), 
3x 

u(L,  t)  =  0, 
du 

u(x,  0)  =  —(x,  0)  =  0. 


forx,-_i  <  x  <  x,  i  =  1, . . . ,  m,  and  f  >  0, 

2r  2r 

for  (n  —  1) —  <  t  <  n  — ,  n  >  1,  and  t  >  0, 
m  m 


(1) 


The  functions  u(x,  t)  and  a(x,  t)  represent  the  displacement  and  stress  respectively.  The  layer  interface  between  the  ith  and 
the  ( i  +  l)th  layer  is  located  at  x,-  for  i  =  1 . . . . ,  m  —  1 ;  x0  =  0  and  xm  =  L,  where  L  represents  the  length  of  the  strip.  By 
replacing  the  spatial  variable  x  with  the  travel-time  variable  defined  in  [14,15]  and  below,1 


t  =  £(x)  = 


ds 

c(s)’ 


it  follows  that  j-  =  c(x)  and  our  initial-boundary  value  problem  (1)  becomes 


(2) 


d2u($,  t)  =  3 2u(g,  t) 

3 12  3£2  ’ 

El  du  -  du 

<7(0,  t)  =  —  —  (0,  t)  =  Ej  — (0,  t)  =/(n), 

ci  3?  3£ 

u(t,  t)  =  0, 
du 

u(t.0)  =  — (£,0)  =  0. 


(i  —  1)t  it 

for -  <  ^  <  — ,  i  =  1, . . . ,  m,  and  t  >  0. 

m  m 

2r  2r 

for  (n  —  1) —  <t  <n  — ,  n>  1,  and  t  >  0, 
m  m 


(3) 


Due  to  the  transformation  of  variables  (2),  problem  (1)  is  simplified  to  a  Goupillaud-type  strip  with  length  equal  to  the 
transit  time  r  through  the  strip,  layer  lengths  equal  to  wave  travel  time  equal  to  xm  for  each  layer  in  either  direction, 
and  wave  speed  of  unity  in  each  layer,  as  seen  from  the  first  equation  in  (3).  The  material  properties  of  the  ith  layer  are 
represented  by  E,-  and  o,-  for  the  physical  case  (1)  and  p,  and  E,  for  the  simplified  case  (3).  Having  p,-  =  %  =  -  =  jEipu 

it  follows  that  the  impedance  ratios  a,-  =  =  3-  and  the  stress  values  £i  =  Eiff, 

i  =  1,  2, . . . ,  m  —  1,  remain  unchanged  as  a  result  of  the  transformation  of  variables  (2).  Therefore,  in  all  our  derivations 


i 


Here,  c  =  c(s)  is  the  piecewise  constant  wave  speed  function  taking  a  constant  value  in  each  layer. 
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of  the  stress  solutions,  we  focus  on  the  simplified  case  (3).  As  shown  in  Fig.  1,  the  original  coordinates  of  the  boundary  and 
layer  interfaces  x0  =  0  <  Xi  <  x2  •  •  •  <  x,  <  •  •  •  <  xm  =  I  are  replaced  by  the  new  coordinates  £0  =  0  <  ^  < 

§2  =  7^  <  •  •  •  <  I,-  =  ^  <  ■  ■  -  $m  =  r  for  i  =  1,2 , . . . ,  m.  The  time  variable  t  is  represented  on  the  vertical  axis  and  Tm 
represents  the  (equal)  wave  travel  time  for  each  layer  of  the  m-layered  strip  in  either  direction.  The  inner  vertical  solid  lines 
represent  the  layer  interfaces. 

Due  to  the  continuity  of  stress  and  displacement  at  each  layer  interface,  the  stress  values  are  constant  in  each  square 
and  half-diamond  subregion.  The  unknown  stress  values  are  represented  by  s,(n)  for  i  =  1,2 , ,m  and  n  >  1.  In  a  given 
layer  i,  the  stress  values  alternate  (interlace)  between  s,_i  (n)  and  s,-(n)  for  i  =  1.2,...,  m,  s0(n)  =  f(n),  and  n  =  1,2,.... 
From  the  superposition  of  the  transmitted  and  reflected  waves  at  the  layer  interfaces  and  boundary,  a  system  of  coupled 
first  order  difference  equations  for  the  stress  terms  was  developed  in  [20]  for  the  special  case  of  a  Heaviside  step  in  stress 
loading/(n)  =  p  =  constant  and  n  >  0.  The  present  work  generalizes  this  system  by  replacing  the  Heaviside  step  in  stress 
loading  with  an  arbitrary  discrete  loading  function/(n),  n  >  0,  as  shown  below 


St(n  +  1)  =  -s  i  (n)  + 
Si(n  +  1)  =  — s,(n)  + 


2a  i 


-s2(n)  + 


~f(n+  1), 


1  +  «!  1  +  «1 
2 «;  2 

si+1(n)  +  — — s,_  i  ( n  +  1), 


1  +  at  1  +  a 

sm{n+  1)  =  -sm(n)  +  2sm_i(n+  1), 


for  i  =  2, . . . ,  m  —  1, 


(4) 


subject  to  the  zero-stress  initial  conditions  s,(0)  =  0  for  1  <  i  <  m.  Here  n  >  0  is  a  suitable  time  index  and 
s,(n)  represent  the  stress  terms  for  1  <  i  <  m  (see  Fig.  1).  The  discrete  function /(n),  n  >  0,  represents  the  stress 
loading  applied  at  £  =  0.  Notice  that  by  applying  the  limiting  condition  a,  —>■  0  to  the  recursive  relation  s,(n  +  1)  = 
— Sj(n)  +  -j^:S1+1(n)  +  j^-Sj-i(n  +  1)  for  i  =  m,  we  recover  the  last  equation  in  (4),  corresponding  to  the  fixed  end.  A 
similar  system  of  difference  equations  may  be  obtained  for  other  initial  and  boundary  conditions. 

System  (4)  can  then  be  solved  using  z-transform  methods.2  As  a  result,  the  z-transform  of  (4)  is  expressed  as 


zS:(z)  =  -S,(z)  + 


2a  i 


- S2(z )  + 


z(F(z)-/(0», 


1  +  oq  1  +  oq 

zS,-(z)  =  -Si(z)  +  2a'  SI+,(z)  +  2  zS,_i(z) 

1  +  or,  1  +  O'; 

zSm(z)  =  —Sm(z)  +  2zSm_i(z). 


for  i  =  2, . . . ,  m  —  1, 


After  reorganizing  the  terms  in  (5),  the  linear  system  is  written  in  matrix-vector  form  as 
AmXm  =  bm, 

where  Am  is  a  tridiagonal  matrix  given  as 


An,  — 


(5) 


(6) 


while 


Xm  — 


r)i  = 


z+  1 

—r)iai 

0 

0 

0 

0 

0 

-T]2Z 

z+  1 

-rj2a2 

0 

0 

0 

0 

0 

-m  z 

z  +  1 

-i)3a3 

0 

0 

0 

0 

0 

0 

0  - 

z  +  1 

0 

0 

0 

0 

0 

Z  +  1 

rs,(z)n 

~m z(F(z)  - 

mr 

S2(z) 

0 

’ 

bm  = 

’ 

_Sm(z)_ 

mx  1 

0 

mx  1 

(7) 


fori  =  1. 


.  m  -1,  r)m  =  2. 


1  +  a,- 

Due  to  the  sparseness  of  bm,  the  solution  of  the  linear  system  (6)  is 

'(-l),+1|A1.ir 
(-1)1+2|A,,2| 


Xm  —  Am  bm  — 


r)\z(F(z)  -/( 0)) 

I  Ami 


(— l)1+m|A1,m| 


(8) 


^  The  definition  for  the  z-transform  ofg(n),  n  >  0,  isZ(g(n ))  =  G(z)  =  n  ^or  lzl  >  ^  m  the  complex  plane  (see  [21]). 
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Here  |Am|  and  |  A-j  j|  are  the  determinants  of  Am  and  its  minor  A^i,  for  i  =  1 ,  ,m.  The  following  recursive  relation  holds: 

|Am|  =  (z  +  l)|Am_i|  —  ^i^2aiz|Am_2|,  (9) 

where  Am  corresponds  to  the  m-layered  strip  with  impedance  ratios  «i ,  a2, . . . ,  am- 1 ;  Am_i  corresponds  to  the  remaining 
(m  —  l)-layered  strip  with  impedance  ratios  a2,  oi3, . . . ,  am_i,  obtained  after  the  removal  of  the  first  layer  from  the  m- 
layered  strip;  Am_2  corresponds  to  the  remaining  (m  —  2)-layered  strip  with  impedance  ratios  a3, . . . ,  am_i,  obtained  after 
the  removal  of  the  first  two  layers  from  the  m-layered  strip.  One  can  derive  by  induction  that  the  determinant  |Am|  is  a 
palindromic  polynomial  with  real  coefficients;  i.e.  the  coefficients  in  front  of  zm_J  and  zJ  are  real  and  equal  to  each-other 

for  j  =  0, _ m  and  m  >  1.  The  complex  roots  of  a  polynomial  with  real  coefficients  occur  in  conjugate  pairs,  and  since 

the  coefficients  are  palindromic,  then  the  roots  also  occur  in  inverse  pairs.  Since  the  inverse  pairing  is  not  necessarily  the 
same  as  the  conjugate  pairing,  if  a  root’s  inverse  is  not  its  complex  conjugate,  then  it  must  be  one  of  a  set  of  four  related 
roots  that  satisfy  a  palindromic  quartic  with  real  coefficients.  As  a  result,  each  palindromic  polynomial  of  even  degree  is 
expected  to  factor  into  quadratic  and  quartic  palindromes  with  real  coefficients.  In  Section  2.2  below,  the  quadratic  and 
quartic  factors  of  |Am|  for  1  <  m  <  5,  are  shown  to  satisfy  the  necessary  and  sufficient  conditions  for  all  their  roots  to  be 
distinct  (complex  conjugate)  on  the  unit  circle.  As  a  result,  the  determinant  |  Am  |  can  be  factored  into  quadratic  palindromes 
with  real  coefficients  for  1  <  m  <  5,  as  presented  later  in  (15). 


2.2.  On  the  zeros  of  the  determinant  |Am| 

The  expressions  for  the  determinant  |Am|  using  the  recursive  relation  (9)  for  1  <  m  <  5  are: 


|Ai  |  —  z  +  1. 
|A3|  =  (z  +  1) 
|As|  =  (z+l) 


,2  _  2  2(Xi-2) 


|A2|  =  (z  +  1 Y  -  =  z- 


Xi 


z  +  1, 


z2-^^z+l 

X2 


4/1 


3,3 


|A4|  =  z - -zJ  + 

X3 


2(4r3  -  X3  +  8) 

X3 


2  4/3 

z2 - -z  +  1, 

X3 


(10) 


4 A  3  2(4_T4  —  X4  +  8)  2 


z’ - z  + 

X4 


X4 


z2 - -z  +  1 

X4 


The  design  parameters  involved  are:  Xm-i  =  nil  0  +  <+)  for  m  >  2,  r3  =  ct\a3  —  1,  and  r4  =  a\a3a4  +  a^a2a4  + 
a\a4  +  a2a4  +  a ^a3  —  1.  Next,  the  necessary  and  sufficient  conditions  (I)  and  (II)  are  applied  to  the  quadratic  and  quartic 
factors  of  |Am|  respectively,  to  prove  that  all  the  zeros  of  |Am|  in  (10)  are  on  the  unit  circle  and  distinct  for  1  <  m  <  5.  As 
for  the  linear  factor  (z  +  1)  which  appears  for  m-odd,  the  zero  z  =  —  1  already  is  on  the  unit  circle  and  distinct  from  all  the 
other  zeros. 

(I)  The  zeros  of  a  quadratic  factor  with  real  coefficients  of  the  form  z2  +  /xz  +  1  are  distinct  ( complex  conjugate)  on  the  unit 
circle  iff  /x2  —  4  <  0  or  equivalently  iff  |/x|  <  2. 

For  m  =  2,  3,  we  have  from  (10)  that  fx  =  (xm  =  —  1(Xm~'~2\  Therefore  the  condition  |/r|  <  2  becomes  equivalent  to 


Xm—  1  3 


<  1,  which  is  true  because  /„ 


>  1. 


(II)  The  zeros  of  a  quartic  factor  with  real  coefficients  of  the  form  z4  +  /xz 3  +  vz 2  +  fxz  +  1  are  distinct  (complex  conjugate) 
on  the  unit  circle3  iff  [/x2  —  4v  +  8  >  0  and  \  |  <  2]. 


For  m  =  4,  5,  we  have  from  (10)  that  /x  =  [xn 


and  v  =  vn 


2(4rm_1-xm-l+8) 

Xm— 1 


The  requirement  fx2  — 4v  +  8  >  0  becomes  equivalent  to  [(xm-i  —  fm-1)2  —  4Xm-i]  =  [(Tm_i  —  2)2  —  4(l  +  rm_i)j  >  0, 
where  Tm_  1  =  Xm-i  —  rm-i- The  desired  inequality  follows  from  the  following  relations: 

(T3  —  2)2  —  4(1  +  r3)  >  (ai  +  a3)2  —  4ai«3  =  (ai  —  a3)2  >  0  for  m  =  4,  and  (T4  —  2)2  —  4(1  +  r4)  > 
(oil  +  q!2  +  ci\Ot2  +  a3  +  a4  +  a3o;4)2  —  4(1  +  r4)  >  (—a^  —  a2  —  oqa^  +  a3  +  a4  +  a3o;4)2  >  0  for  m  =  5. 


The  requirement 


— n  it  4y+8 

2 


<  2  is  equivalent  to 


On— l=t^/  (Xm— 1— On— 1)^— 4Xm— 1 

Xm— 1 


<  1,  which  is  true  because 


Xm—  1  >  1  >  0  and  rm_  1  +  1  >  0  for  m  =  4,  5.  As  seen  later  in  (A.5)-(A.7),  the  two  requirements  in  (II)  guarantee 
real  cosine  values. 

For  the  general  m-layer  case,  our  preliminary  results  with  tridiagonal  Toeplitz  matrices,  show  that  we  can  find  designs 
with  any  number  of  layers  for  which  all  the  zeros  of  |Am|  are  distinct,  have  no  multiples,  and  are  on  the  unit  circle.  Indeed, 


3  This  criterion  is  derived  after  dividing  the  quartic  polynomial  by  z2  and  then  applying  the  transformation  y  =  z  +  - . 
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using  standard  techniques,  the  symmetric  matrix  Gm  can  be  formed  such  that  |Am|  —  |Gm|, 

•••0  0  0  0 


"z  +  1  Vz£  1  0 

yfzb  z+  1  yjzi; 2  0 

0  yj  Zt,2  Z  +  1  yfzX 3  0 


0  vjz^m— 2  Z  1 

0  0  y[zU^\  z+  1  . 


(11) 


with  £  =  aii)i?]i+i  for  i  =  1,  2, . . . ,  m  —  1.  As  before,  ??,  =  for  i  =  1, . . , ,  m  —  1  and  =  2.  If  Gm  has  constant 

off-diagonal  terms,  i.e.,  VzfT  =  =  ■  ■  ■  -Jzbflf  =  ~Jzt,  or  equivalently  t,  =  for  i  =  1,  2, . . . ,  m  —  1.  Then  Gm 

is  a  tridiagonal  Toeplitz  matrix  with  palindromic  determinant  |Gm|  which  can  be  written  as  [22], 


m 

|Gml  =  n  (z  +  1  “  m))  , 


(12) 


k=l 

or  equivalently, 


IGml 


1 1  [z2  +  2z  (l  —  2 fg2(k,  m))  +  l]  for  m  even, 

k= 1 

L?J 

(Z+l)  n  [z2  +  2z  (l  —  2fg2(k,  m))  +  l]  form  odd, 


(13) 


where  g(k,  m )  =  cos  (^j).  For  0  <  K  <  1,  the  zeros  of  the  quadratic  factors  in  (13), 


zk  =  2fg2(fc,  m)  -  1  ±  (14) 

are  distinct,  complex  conjugate,  and  on  the  unit  circle,  with  |z/(|2  =  ZyZ\<_  =  1  and  arglz^)  =  0/(,  for  all  1  <  k  <  [f  J-  For 
these  f-values,  the  representation  in  (13)  becomes  consistent  with  (15).  In  our  numerical  experiments  displayed  later  in 
Fig.  5a,  b,  for  a  given  value  of  f ,  0  <  f  <  1,  the  impedance  ratios  are  derived  from  the  recursive  relations:  am_i  =  and 
a,  =  —  1  +  4_(J(a.+i  for  i  =  1, . . . ,  m  —  2.  Other  tridiagonal  [23,24]  or  k-Toeplitz  [25,26]  matrices,  which  have  analytical 
expressions  for  their  determinants  like  that  given  by  ( 12),  may  also  be  studied  in  an  attempt  to  generalize  the  class  of  media 
for  which  the  roots  of  |Gm|  =  0  are  distinct  and  on  the  unit  circle  for  arbitrary  m. 

In  this  paper,  we  only  consider  m-layered  designs  for  which  all  the  zeros  of  |Am|  are  distinct,  have  no  multiples,  and  are 
on  the  unit  circle.  Based  on  the  results  above,  this  class  of  designs  is  not  empty.  For  z^  =  el0k  =  cos  9^  +  1  sin  0/(,  we  have 
that  z,7 1  =  fit  =  e10k  =  cos  9^  —  I  sin0/(  and  (z  —  e1Bk)(z  —  e~16k)  =  z2  —  2z  cos  9k  +  1.  As  a  result,  a  newly  factored 
representation  of  |Am|  is  obtained,  with  palindromic  quadratics  with  real  coefficients  involving  the  cosines  of  the  angles 
9\,  92, . . . ,  9^m j  and  linear  factor  (z  +  1)  for  m-odd, 


I  Am  I  = 


L?J 

[z2  —  2z  cos  9k  +  1] 

k=l 

LfJ 

(z  +  1)  PJ[z2  -  2 z  cos  9k  +  1] 

fc=i 


for  m  even, 


for  m  odd. 


(15) 


2.3.  General  formulas  for  the  stress  and  resonance  frequencies 

Based  on  the  factorization  of  the  determinant  |Am|  given  in  (15),  for  any  m-layered  design,  if  there  are  m  distinct  roots 
on  the  unit  circle,  there  are  m  distinct  angles.  The  [y]  essential  angles  0  <  9k  <  n  for  k  =  1, . . . ,  LyJ  were  known  in  [20] 

i.j 

as  the  base  angles.  The  equations  for  the  angles  90  =  jz  and  {9k}kl\  in  terms  of  the  design  parameters  (combinations  of  the 
impedance  ratios)  for  up  to  five  layers  are  derived  in  [20],  Since  the  degree  of  |Am|  is  m  and  the  degree  of  |Ai  j|  is  m  —  1, 
for  i  =  1,2 , . . . ,  m,  the  substitution  of  (15)  into  (8),  results  in  the  following  expansion  of  the  components  xm(i)  ofxm  into 
partial  fractions: 
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Xm(0  =  Z(F(Z)-f(  0)) 


=  z(F(z)-/(  0)) 


4i(-1)'|Aij| 


I  Am  I 


LfJ 


bi.o _ a*kz  +  b*k 

1  tlz2~2: 


z  4- 


2  z  cos  0k  +  1 


=  (F(z)  —/CO)) 


I  —  I 

bi.oz  o*  kz2  +  b*„z 


z  + 


:,(K  y- _ 

-fl  z2  - 


2z  cos  0*  +  1 


(16) 


where  bi0  =  0  for  m  even  and  i  =  1,2, ...  ,m.  After  applying  the  inverse  z-transform  to  (16),  a  general  representation  for 
the  stress  terms  in  (4)  is  obtained, 


Si(n)  =  Z  1  (xm(i))  =  /(n)  * 


LfJ 


bi,o(-l)n  +  E  ai,k  cos(n&i<)  +  bUk  sin (n9k) 


k=  1 


— /(0)  • 


M-1)"  +  E  aukcos(n6k)  +  bLk  sin(n9k) 


(17) 


for  the  following  choices  of  the  coefficients, 
a*kcos9k  +  b*k 


@i,k  —  ^i,k’ 


hk  = 


for  k  = 


- . m 


sin  0k 

for  i  =  1,2, ... ,  m  and  n  >  0.  The  operation  *  means  convolution.  The  stress  representation  previously  derived  in  [20] 
can  now  be  recovered  for  the  special  choice  of  loading/(n)  =  p  for  n  >  0.  The  stress  representation  in  (17)  involves  sums 
of  sines  and  cosines  of  multiples  of  angles,  which  relate  to  the  Chebyshev  polynomials  of  the  first  and  second  kind.  For  the 
choice  of f(n)  =  sin(na>),  n  >  0,  after  substituting /(0)  =  0  and  F(z)  =  Z(f(n))  =  Z( sin(n<y))  =  (z2 _z2zcos>1s>+ 1 1  'nt0  ^6), 
and  then  applying  the  inverse  z-transform  we  get 


s,(n)  =  bt0Z  1 


z2  sin(ftj) 


(z2  —  2z  cos  a>+  1 )  (z  +  1) 


LfJ 

+  Ez-' 


z  sin(w)  •  ( a*kz 2  +  b*kz) 


(z2  —  2z  cos  9k  +  1  )(z2  —  2  z  cos  w  +  1) 


(18) 


In  the  subsections  that  follow,  the  inverse  z-transform  is  applied  to  the  expressions  in  (18)  to  derive  the  stress  formulas 
at  resonance  and  non-resonance  frequencies;  the  results  are  verified  with  numerical  experiments.  Based  on  these  results, 
we  prove  that  the  resonance  frequency  spectrum  for  the  discrete  model  consists  of  all  the  positive  angles  coterminal  with 

L®J 

0O  =  n  or  {±0/;}^! ,  as  shown  below 


wi0  =  0o  +  2/07r  =  (2 10  +  1)jt  (for  m  odd  only), 

=  +  2/iJ t ,  lo,  l\,  h  =  0.  1.  2,  . . 

&>i2,k  —  —  9k  +  2(12  +  1)jt,  k  =  1, . . . ,  |_yj 
or  equivalently, 


cosw  : 

COS  ft)  : 


—  1  (for  m  odd  only), 
cos  0k  for  k  =  1, . . . ,  |^— J  . 


(19) 


(20) 


The  fact  that  the  0-angles  depend  only  on  (certain  combinations  of)  the  impedance  ratios  implies  that  the  resonance 
frequencies  for  the  discrete  model  in  (19)  depend  on  the  same  combinations  of  the  impedance  ratios  and  no  other 
parameters. 


2.4.  Stress  solutions  at  non-resonance  frequencies 

Here  we  show  that  the  values  of  the  driving  frequency  w  that  satisfy  cos  w  ^  cos  0O  =  —  1  (for  m  odd),  and  cos  w  ^  cos  9k 
for  all  k  =  1, . . . ,  LfJ,  represent  non-resonance  frequency  values.  In  the  set  of  real  numbers,  this  set  of  frequency  values  is 
the  complement  of  the  set  of  values  identified  in  (20).  Indeed,  applying  the  partial  fraction  expansion  to  (18),  we  derive  the 
following  stress  solutions, 


bj  o  sin  cw 

=  -477^ - ^r(-l)n  + 


2(1  +  cos  &)) 


bko  sin  a) 
2(1  +  cos®) 


LtJ 

+  E At 


cos(nftj) 


+ 


b  L5J 


LfJ  LfJ 

sin(nftj)  +  y  Q,k  cos (n9k)  +  ^  Di<k  sin(n0fc), 


(21) 
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Fig.  2.  Stress  time  history  for  a  two-layered  Goupillaud-type  strip  in  the  middle  of  the  second  layer  located  at  £  =  3/4.  Impedance  ratio  a i  =  1  /3,  angle 
0 !  =  2tt /3,  t  =  1,  and  loading/(n)  =  sin(na))  at  £  =  0.  (a)  &>  =  7r /4.  (b)  co  =  0.9  •  =  0.67T. 


where 


=  ~Q,k  = 


2(cos  w  —  cos  Ok) 


a*kcosw  +  b*k 
2(cos  «  —  cos  8k) ' 


(a*k  cos  8k  +  b*k)  sin  a) 
2  sin  0/((cos  a;  —  cos 


(22) 


As  before,  bi0  =  0  when  the  number  of  layers  m  is  even.  For  a  given  value  of  w,  the  stress  solutions  in  (21),  including  all 
the  coefficients  and  angles  8k  for  k  =  1, . . . ,  [  ,  can  be  also  expressed  in  terms  of  the  impedance  ratios  only  (see  (6)-(8)). 

Illustrations  of  the  bounded  stress  solutions  at  non-resonance  frequencies  appear  in  Fig.  2.  As  depicted  in  the  plots  in  Fig.  2, 
when  positioned  in  the  middle  of  the  second  layer,  the  time  interval  at  which  s,  (n)  is  reached  is  ^ r ~  <  r  ±  for 

n  >  1  and  i  =  1,2  (see  Fig.  1).  Beat  phenomena  are  clearly  visible  in  2b,  as  the  driving  frequency  approaches  a  resonance 
frequency. 


2.5.  Stress  solutions  at  resonance  frequencies 


The  stress  solutions  in  (21)  suggest  that  the  driving  frequencies  co  identified  in  (20)  are  resonance  frequencies.  Flere  we 
confirm  this  by  deriving  the  stress  solutions  at  these  expected  resonance  frequencies,  demonstrating  that  such  solutions 
grow  without  bound  over  time  while  oscillating. 

Case  I.  cos(ftj)  =  —  1  =  cos  80  (for  m  odd  only). 

The  fact  that  cos(<J>)  =  —1  or  equivalently  that  a>  takes  values  cbi0  =  (2 10  +  l)7r  for  l0  =  0,  1,  2, . . .,  implies  that 
/ (n)  =  sin(na))  =  sin(n7r)  =  0  for  all  n  >  0.  As  a  result,  due  to  the  zero  input  in  loading,  we  will  not  be  able  to  detect 
resonance,  as  concluded  by  (17).  By  choosing/(n)  =  cos  (nob)  =  cos(njr)  we  overcome  this  obstacle.  Based  on  (16)-(17), 
we  expect  any  unbounded  stress  terms  to  come  from 


F(z)- 


bj.oz 
z  + 


m  +  Kkz 

z 2  —  2z  cos  8k  +  1 


T  +  S 


Z 

z  +  1 


bi,0z 


LtJ 


t+£ 


“tkz2  +  K,<z 


z  +  1  z2  —  2z  cos  &k+  1 


(23) 


after  the  inverse  z-transform  is  applied.  We  conclude  that  the  only  term  in  (23)  that  generates  unbounded  stress  values  is 


bj,  QZ 

(z+ty 


for  bj  o  0,  due  to  the  fact  that 


Z-1 


bj.oz2 
(z  +  l)2 


=  £>i,o  •  VV„(— 1)  =  £),-  0(— l)"(n  +  1). 


(24) 


Here  Wn(y)  represents  the  Chebyshev  polynomial  of  the  second  kind  evaluated  at  y  =  —1.  Notice  that  the  resonance 
frequency  spectrum  =  (2 10  +  1)jt  for  l0  =  0,  1,  2, . . .,  is  universal  for  all  the  designs  with  an  odd  number  of  layers 
and  it  is  independent  of  any  design  parameters.  This  is  illustrated  in  Fig.  3  for  the  seven-  and  eleven-layer  case. 

Case  II. a.  cos  a>  =  cos  8k*  and  sin  to  =  sin  8k*. 
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Fig.  3.  Stress  time  history  for  a  multilayered  strip  with  r  =  1  and  impedance  ratios  a\  =  3,  a2  =  2,  a3  =  1.5,  a4  =  2.2,  a5  =  0.3,  a6  =  1.7,  a-]  =  1.4, 
as  =  3.1,  ofg  =  0.8,  ofio  =  4  as  applicable;  loading/(n)  =  cos(no>)  is  applied  at  £  =  0  for  co  =  n.  (a)  The  middle  of  the  second  layer  of  a  seven-layered 
strip  located  at  f  =  3/14.  (b)  The  middle  of  the  second  layer  of  an  eleven-layered  strip  located  at  f  =  3/22. 


Suitable  partial  fraction  expansions  of  ( 1 8 )  give  the  following  stress  solutions  for  the  frequency  spectrum  cbixk*  =  9k*  +  21]  tt  , 

Zi  =0.  1.2,..., 


si(n) 


bj  o  sin  9k* 

— - - - — (-1)"  + 

2(1  +  cos  Ok*) 


b,- o  sin0/(* 
2(1  +  cos  6k*) 


Lf  J 

+  GlAik*  +  A  ,k 

k=\,k^k* 


cos  (ndk*) 


+ 


bi,0 


L?J 


— - 1"  nBi,k*  +  Q 'k*  + 

^  k=l,tyk* 


sin  (ndk*) 


Lf J  Lf J 

+  Y  C„fc  cos(nOk)  +  Y  DLksm(n9k).  (25) 

(t=l  ,tyk*  k=l,tyk* 

Here  Bi  k,  Qik,  are  given  from  (22)  for  w  =  6k *  +  2/jjr,  /j  =  0,  1,  2, . . .,  while 


A,k*  =  — 


(a*k*  cos  9k,  +  b*l(t) 
2  sin  9k* 


Bu 


i,k* 

~Y' 


_  aik*  +  Kk *  C0S  °k* 
2  sin2  9k* 


(26) 


As  before,  bj,0  =  0  for  m  even  and  sinBj-  ^  0  for  all  k  =  1, . . . ,  [  jJ .  All  the  terms  in  the  stress  solution  (25)  are  bounded, 
except  cos(n9k *)  +  Bj k*  sin(n0/<*)]  •  n,  which  are  multiplied  by  the  time-related  linear  factor  n.  This  implies  that  the 
stress  amplitude  is  expected  to  grow  without  bound  while  oscillating  over  time,  proving  resonance. 

Case  Il.b.  cos  cy  =  cos9k*  and  sin«  =  —  sin 0^*  =  sin(— 9k»). 


A  similar  stress  representation  to  (25)  can  be  obtained  at  the  frequency  values  a>i2tk*  =  —0k*  +  2 (l2  +  1)tt ,  l2  =  0,  1,  2 . 

after  a  sign  adjustment  for  some  of  the  coefficients. 

The  resonance  response  in  a  four-  and  five-layered  elastic  strip  using  the  analytical  results  generated  from  (25)  is 
illustrated  in  Fig.  4.  This  is  a  demonstration  and  corroboration  of  formula  (19)  in  predicting  resonance.  As  illustrated  in 


Fig.  4  for  the  four-layer  case,  the  base  angle  92  =  cos  1 


C3-V  (X3-C3)2-4/3 

X3 


represents  a  resonance  frequency.  The  same 


is  true  for  the  five-layer  case  and  the  base  angle  9\  =  cos  1 


|  C4+x/(/4-/4)2-4x4 

V  X4 


.  The  relevant  design  parameters  are  given 


by  Xm- 1  =  ITwi1^  +  ai)  f°r  m  =  4,  5;  r3  =  oqa3  —  1,  and  r4  =  oqaqa^  +  oqa 2aA  +  +  a2a4  +  aqo!3  —  1  while  the 

impedance  ratios  are  =  0.6,  a2  =  1.5,  a3  =  2,  and  a4  =  1.2,  as  applicable. 

The  resonance  response  in  a  twenty  two-  and  thirty  nine-layered  elastic  strip  generated  from  the  recursive  relations  (4) 
and  (13)-(14)  is  illustrated  in  Fig.  5. 

This  is  another  demonstration  and  corroboration  of  formula  (19)  in  predicting  resonance.  For  the  twenty  two  layer  case, 
the  base  angle  corresponding  to  k  =  11  and  f  =  7/9  represents  a  resonance  frequency.  Similarly  for  the  thirty  nine  layer 
case,  the  base  angle  obtained  for  k  =  5  and  £  =  1/2  represents  a  resonance  frequency.  The  impedance  ratios  in  both  cases 
are  derived  from  the  recursive  relations:  am_i  =  -p—  and  cq  =  —  1  +  .  4 - for  i  =  1, . . . ,  m  —  2. 

4  5  4  5  5ai+ 1 
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Fig.  4.  Stress  time  history  for  a  multilayered  strip  with  r  =  1  and  loading /(n)  =  cos  (neb)  applied  at  f  =  0.  (a)  The  middle  of  the  second  layer  of  a 
four-layered  strip  located  at  f  =  3/8,  when  the  driving  frequency  is  cb  =  02, 02  =  cos-1  ^r'3~V/(x3-^3)  ~4*3  ^  ^  _  12,  r3  =  0.2.  (b)  The  middle  of  the 

second  layer  of  a  five-layered  strip  located  at  £  =  3/10,  when  the  driving  frequency  is  cb  =  0\,  0\  =  cos-1  (  Fa+^(u  r<4>  4x4  \  =  26.4,  F4  =  5.24. 


Fig.  5.  Stress  time  history  for  a  multilayered  strip  with  r  =  1  and  loading/(n)  =  sin(n£>)  applied  at  f  =  0.  (a)  The  middle  of  the  first  layer  of  a  twenty 
two-layered  strip  when  the  driving  frequency  is  cb  =  On,  (b)  the  middle  of  the  fifth  layer  of  a  thirty  nine-layered  strip  when  the  driving  frequency  is 

cb  =  6$. 


2.6.  Accuracy  and  precision  of  the  recursive  and  explicit  solutions 


Despite  the  fact  that  the  recursive  relations  (4)  and  the  explicit  formulas  (17)  are  exact,  the  stress  solutions  s,(n),  n  >  1, 
i  =  1 , ,m,  calculated  from  these  formulas,  are  not  always  exact.  Here  we  show  that  for  a  Heaviside  stress  loading  and 
certain  layer  properties,  the  stress  solutions  represented  by  terminating  decimals  are  exact,  as  they  do  not  exhibit  roundoff 
error.  In  others  cases,  when  the  stress  solutions  are  represented  by  non-terminating  decimals,  a  precision  loss  is  observed. 
This  discussion  is  particularly  important  for  verification  of  numerical  finite  element  solutions  with  “exact"  solutions;  one 
should  be  aware  of  the  loss  in  precision  even  in  the  “exact”  solutions  [27], 

For  a  free-fixed  two-layered  strip  subjected  to  a  Heaviside  loading  with  layer  impedance  ratio  a  =  3,  the  stress 
in  layer  2  is  represented  by  a  periodic  sequence  of  12  rational  numbers  that  can  be  expressed  as  terminating  decimals 
{ ’ ,  1,  |,  2,  2,  2,  1,  0,  0,  0, . . .}  (see  Fig.  7(b)).  Both,  the  recursive  relations  (4)  and  the  explicit  trigonometric  formulas 

(17)  predict  this  stress  sequence  exactly  for  all  n  >  1. 

Alternatively,  for  a  free-fixed  two-layered  strip  subjected  to  a  Heaviside  loading  with  layer  impedance  ratio  a  =  4,  a  loss 
in  precision  is  observed  when  a  fixed  number  of  significant  digits  is  used  in  the  numerical  calculations.  The  loss  of  precision 
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Fig.  6.  Precision  loss  in  stress  calculations  for  a  two-layered  Goupillaud-type  strip  with  a  =  4,  subjected  to  a  Heaviside  loading /(n)  =  p.  (a)  Linear  loss 
in  precision  obtained  from  the  recursive  relations  (4),  and  (b)  log-linear  loss  in  precision  in  the  explicit  formula  (17). 


with  the  recursive  relations  (4)  is  due  to  roundoff  error,  while  the  loss  of  precision  with  the  trigonometric  formulas  (17)  is 
due  to  a  loss  of  significant  digits  in  the  trigonometric  evaluations  that  involve  increasingly  large  arguments. 

Using  arbitrary  precision  arithmetic  available  in  the  commercial  software  package  Mathematica  8,  the  loss  of  precision 
in  stress  calculations  is  illustrated  in  Fig.  6.  If  the  calculations  begin  with  numerical  precision  of  500  decimal  places,  the 
recursive  relations  exhibit  a  linear  loss  in  precision  with  each  recursive  step,  whereas  the  explicit  analytical  solutions 
exhibit  a  log-linear  loss  of  precision.  Similar  conclusions  can  be  made  for  the  recursive  and  analytical  stress  solutions  at 
the  resonance  and  non-resonance  frequency  cases  studied  here. 


3.  Continuous  forcing  function  and  natural  frequencies 


In  this  section,  we  consider  the  initial-boundary  value  problem  (3)  for  a  stress  loading  of  the  form  cr  (0,  t)  =  sin(a>t),  for 
f  >  0.  Unlike  the  discrete  model,  the  stress  loading  (forcing  function)  for  the  continuous  model  is  continuous  with  time.  In 
addition,  unlike  the  discrete  model,  the  resonance  frequency  spectrum  for  the  continuous  model  is  not  found  by  solving  the 
corresponding  initial-boundary  value  problem.  Instead,  we  use  the  linear  transformation  in  frequency 


w  =  — a>,  (27) 

m 

between  the  frequency  of  a  continuous  forcing  function  a>  and  the  frequency  of  the  corresponding  discrete  forcing  function 
a>  and  time  step  to  extend  the  resonance  frequency  results  (19)  to  our  continuous  model: 


«;0  = 


m 

2r 

m 

2t 

m 

2t 


Wo  +  2/q7T] 


1 


(2/q  +  1  )tt  , 


H\) 

Wk  +  2/iJr]  =  ,  .  •  Wk  +  2/i 7r ] , 


2  a) 

•  [~6k  +  2(12  +  1)?f]  = 


1 


2  a) 


•  [~6k  +  2(h  +  ] , 


(for  m  odd  only), 

Iqi  1i>  I2  =  0>  1  >  2, . . . , 


(28) 


The  linear  transformation  in  frequency  (27)  is  discussed  in  [28],  as  the  nature  of  time  (continuous  or  discrete)  is  expected 
to  affect  the  nature  of  the  frequency.  The  frequency  for  the  analog  (continuous)  time  signal  is  related  to  the  frequency  of 
the  digital  (discrete)  time  signal  through  the  linear  transformation  (27),  involving  the  time  step  or  the  so-called  sampling 
interval.  For  our  problem,  one  can  derive  relation  (27)  between  the  frequency  co  of  the  continuous  forcing  function  sin(a>t), 
f  >  0,  and  the  frequency  w  of  the  discrete  forcing  function  sin(am),  n  >  1,  based  on  the  following  operations, 


sin(ta>)  =  sin 


sin 


t  I  2t  \  /  2r  \ 

=-  •  — co  1  =  sin  n  ■  — co  =  sin(ntt>). 

If  m  )  V  rn  ) 


As  seen  in  Fig.  1,  t  >  0  represents  the  time  variable,  n  =  [“tt!  =  0,  1,  2 . represents  a  time-related  index,  and  ^ 

TrT 

represents  the  (equal)  wave  travel  time  for  each  layer. 

The  proposed  resonance  frequency  spectrum  in  (28)  also  represents  the  natural  frequencies  of  a  free-fixed  layered 
strip.  The  proposed  natural  frequency  formulas  in  (28)  suggest  that  for  a  free-fixed  m-layered  Goupillaud-type  strip,  the 
natural  frequency  spectrum  depends  on  (certain  combinations  of)  the  impedance  ratios  a-i,  a2,  ■  •  • .  «m-i  and  it  is  inversely 
proportional  to  the  (equal)  wave  travel  time  T  for  each  layer. 
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Fig.  7.  Stress  time  history  for  a  two-layered  Goupillaud-type  strip  with  impedance  ratio  a \  =  3, 6\  =  | ,  r  =  1,  at  the  middle  of  the  second  layer  located 
at  £  =  3/4  subjected  to  loading  (a )/(n)  =  sin(n<y)  with  co  =  6  =  j.  (b)/(n)  =  1  for  n  >  0. 


3.1.  An  alternative  derivation  of  the  natural  frequency  spectrum  using  the  frequency  equation 


Here,  we  derive  rigorously  the  natural  frequency  spectrum  (28)  for  a  free-fixed  m-layered  Goupillaud-type  strip,  by 
generalizing  the  procedure  described  in  [10]  and  solving  the  frequency  equation  for  1  <  m  <  5.  Similar  to  Section  2,  the 
transformation  of  variables  (2)  is  used  to  replace  the  spatial  variable  x  with  the  travel-time  variable  £  for  the  boundary  value 
problem  under  consideration.  The  advantage  of  applying  the  transformation  of  variables  (2)  is  that  the  frequency  equation 
obtained  is  a  lot  simpler  to  solve,  allowing  us  to  recognize  patterns  and  develop  formulas  for  the  natural  frequency  spectrum. 

Following  the  approach  taken  in  Section  2,  and  applying  the  transformation  of  variables  (2),  the  governing  equations  for 
u(§,  t)  become: 


d2Uj($,t)  =  d2u ,-(g,  t) 
dt2  d%2 


for  i  <  f  <  |i,  i=l,...,m, 


whereto  =  0  <  ^  <  £2  =  7^  <  •••  <  ^  =  r,  for  i  =  1,  2,  ...,  m.  Notice  that  the  wave 

speed  has  become  unity  for  each  layer,  while  the  material  properties  become  p,  =  £,  =  ^  =  y/Efpi  for  i  =  1,2, ....  m, 
as  previously  defined  in  Section  2.  Seeking  harmonic  solutions,  we  may  represent  u(§,  t)  =  U(t()if(t)  and  [/,•(£)  = 
Aj  sin  cot;  +  Bj  cos  co%,  for  |,_i  <  f  <  §,•  and  i  =  1 , ...  ,m,  where  co  is  the  frequency  of  the  harmonic  motion.  Similarly, 
if(t)  can  be  represented  in  the  form  ifr(t)  =  A*  sin  cot  +  B*  cos  cot  or  in  exponential  form  as  \[r  (t)  =  A**ela>t  +  B**e~lcot, 
where  2  =  v7— T.  The  boundary  conditions  at  £  =  0  and  £  =  r  imply  that 


dt/jfO) 

d| 


=  Um(r)  =  0, 


which  together  with  the  continuity  of  stress  and  displacement  conditions  at  each  layer  interface,  lead  to  the  following 
homogeneous  linear  system  of  equations, 


Ax  =  0, 

Aj  sin  (co£,)  +  Bj  cos  (cof,-)  =  Ai+ 1  sin  (cofj)  +  Bi+i  cos  (co£ ;)  for  i  =  1 , ...  ,m  -  1 , 

a;  [4j  cos  (co£i)  —  Bj  sin  (co^)]  =  A+i  cos  (co^)  —  Bi+1  sin  (co^)  for  i  =  1 , ...  ,m—  1, 

sin  (cot;m)  +  Bm  cos  (cu£m)  =  0. 


(29) 


The  relations  ^  =  07,  for  i  =  1, . . . ,  m  —  1,  were  applied  to  the  system  (29),  using  the  fact  that 

c  ■  p  =  E/c.  The  system  (29)  is  written  in  matrix-vector  form  as 

4mv  =  6, 


(30) 


G.A.  Gazonas,  A.P.  Velo  /  Wave  Motion  49  (2012)  135-151 


147 


where 


A  in  — 
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(31) 


and 
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—  COS  Y2 

0 

a2  cos  Y2 

—ol2  sin  Y2 

—  COS  Y2 

sin  Y2 

— 

0 

0 

0 

0 

0 

sin  Ym-i 

COS  Ym—  1 

-  sin  Ym-i 

COS  Ym—  1 

0 

0 

dm- 1  COS  Ym-i 

— ttm-1  sin  Ym-\ 

-  COS  Ym-1 

sin  Ym-i 

0 

0 

0 

0 

sin  Ym 

COS  Ym 

(32) 


with  Yi  =  ^  for  i  =  1, . . . ,  m.  Seeking  a  nonzero  solution  for  the  homogeneous  system  (30),  the  determinant  of  the 

system  matrix  Am  must  vanish 


\Am\  =0. 


(33) 


Eq.  (33)  gives  the  frequency  equation  for  co,  which  can  be  successfully  solved  using  symbolic  algebra  software  for  a  strip  with 
up  to  five  layers  (see  Appendix  A).  Solving  the  frequency  equation  (33)  for  the  general  m-layer  case,  poses  a  computational 
challenge. 

Comparing  the  formulas  for  the  natural  frequencies  from  Appendix  A  with  the  formulas  for  the  0-angles  from  [20],  we 
obtain  the  following  relation, 

2r 

cos  — co  =  cos  0, 
m 

which  confirms  the  relationship  in  (27)  between  the  frequencies  co  given  in  (19)  and  co  in  (28). 


4.  Applications  of  the  frequency  results 

The  properties  of  the  materials  used  in  this  section  are  included  in  Appendix  C.  The  characteristic  impedance  in  each  layer 
is  calculated  in  terms  of  the  material  properties  s/Ep.  Once  the  length  of  one  layer  is  chosen,  the  length  of  each  remaining 

layer  is  determined  according  to  the  recursive  relations:  -  =  —  or  Lj+i  =  / ^±1  •  L,,  to  provide  equal  wave  travel  time. 

ci  ci+ 1  V  Pi+ 1 

Here  Lj,  Cj,  Ej  and  pj,  for  j  =  i,  i  +  1,  represent  the  length,  wave  speed,  elastic  modulus  and  material  density  for  the  ith  and 
(i  +  1 ) th  layer  respectively,  i  =  1, . . . ,  m  —  1. 


4.1.  One  dimensional  layered  media  with  a  common  frequency  spectrum 

According  to  ( 19),  the  resonance  frequency  spectrum  for  the  discrete  model  depends  only  on  certain  combinations  of  the 
impedance  ratios  a\ ,  a2, . . . ,  am_i .  Such  combinations  of  impedance  ratios  when  2  <  m  <  5  are  represented  by  the  design 
parameters  Xm- i  f°r  m  =  2,3,  and  Xm-i.  I'm- i  f°r  m  =  4,  5,  see  (10).  There  are  (theoretically)  infinitely  many  m-layered 
designs  of  a  Goupillaud-type  elastic  strip  that  have  the  same  values  of  the  design  parameters  Xm-i.  I'm- i  for  2  <  m  <  5, 
and  hence  the  same  resonance  frequency  spectrum  for  the  discrete  model. 

As  for  the  natural  frequency  spectrum  of  a  free-fixed  m-layered  Goupillaud-type  strip,  based  on  (28),  it  depends  not 
only  on  combinations  of  the  impedance  ratios  a^,  a2,  •  •  • ,  <*m-i  but  also  on  the  (equal)  wave  travel  time  T  for  each  layer. 
Therefore,  to  obtain  designs  with  a  common  natural  frequency  spectrum,  besides  having  the  same  values  of  the  relevant 
design  parameters,  the  layer  lengths  have  to  be  chosen  so  that  they  all  have  the  same  wave  travel  time. 
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For  instance,  for  the  three-layer  case  with  /2  =  (1  +  oq)(l  +  a2 ),  the  tungsten-lead-aluminum  configuration  with 
oq  ~  5.828  and  a2  ~  0.952  has  the  same  value  of /2  ~  13.33  as  the  lead-aluminum-copper  configuration  with  oq  ~  0.952 
and  (*2  ~  5.828.  Another  material  design  with  X2  ~  13.33  is  the  aluminum-steel-brass  configuration  with  oq  ~  0.331  and 
a2  ~  9.016.  In  addition,  in  order  for  these  material  designs  to  have  a  common  natural  frequency  spectrum  for  the  free-fixed 
strip,  only  the  length  of  one  layer  in  one  of  the  designs  may  be  chosen  at  random.  The  lengths  of  the  remaining  layers  in  all 
the  three  designs  are  then  determined  to  provide  equal  wave  travel  time. 

4.2.  Design  modification  that  gives  a  desired  frequency  spectrum  within  limitations 

Adding  a  new  layer  in  front  of  an  existing  two-layered  Goupillaud-type  strip  with  xi  =  (1  +  oq)  and  9\,  allows  us  to 
modify  the  resonance  frequency  spectrum  for  the  discrete  model  in  (19)  and  (20).  The  length  of  the  new  layer  must  be 
chosen  so  that  the  strip  remains  of  Goupillaud-type.  The  new  three-layered  strip  is  characterized  by  X2  =  (1  +  a*)(l  +  a|) 
and  base  angle  9*,  with  to  be  determined  and  a\  =  oq.  The  unknown  impedance  ratio  a*  is  derived  from  the  following 
formula 

2  xi  —  2 

a.*,  =  —  1  -| - ,  where - =  cosfA  <  V  =  cos  6?  <  1, 

1  Xi(l-V)  Xi  1 

for  a  desirable  value  of  the  cosine  of  9*  given  by  V  =  cos  9*  =  For  instance,  for  a  given  6\  =  j  and  for  the  choice  of 
9*  =  |  we  have  that  cos  6»i  =  cos  (|-)  <  V  =  cos  0*  =  cos  (|),  and  the  impedance  ratios  of  the  three-layered  strip  are 
a*,  =  —  1  +  and  a\  =  The  new  natural  frequency  spectrum  is  obtained  by  substituting  9?  =  f  into  (19) 

or  (20). 

4.3.  Resonance  frequencies  and  optimal  designs 

For  the  m-layered  strip  subjected  to  a  Heaviside  loading  [20],  it  was  shown  that  when  the  optimal  values  of  the  base  angles 
Otopt  for  k  =  1, . . . ,  LfJ,  given  in  Table  B.l  of  Appendix  B,  were  substituted  into  the  respective  angle-design  parameter 
equations,  optimal  m-layered  designs  were  generated  for  2  <  m  <  5.  These  optimal  designs  were  expected  to  provide  the 
smallest  stress  amplitude  of  double  the  loading  in  all  layers,  for  all  time. 

However,  when  these  optimal  angle  values  are  substituted  into  (19)  and  (28),  they  generate  the  resonance  frequency 
spectrum  for  the  discrete  and  continuous  model  respectively.  As  a  result,  for  these  designs,  depending  on  the  type  of  loading 
one  can  get  either  the  best  or  the  worst  outcome  when  it  comes  to  controlling  the  stress  amplitude,  see  the  illustration 
in  Fig.  7.  As  the  stress  amplitude  in  Fig.  7(a)  grows  without  bound  over  time  for  a  harmonic  forcing  condition,  the  stress 
amplitude  in  Fig.  7(b)  does  not  exceed  double  the  loading  for  a  Heaviside  loading  condition. 


4.4.  Natural  frequencies  of  a  free-fixed  non-Goupiliaud-type  layered  strip  with  integer  layer  length  ratios 


The  natural  frequency  results  for  a  free-fixed  Goupillaud-type  layered  strip  with  equal  layer  lengths  given  in  (A.1)-(A,7), 
may  be  extended  to  a  free-fixed  non-Goupillaud-type  layered  strip  with  integer  layer  length  ratios. 

For  instance,  the  frequency  results  for  a  free-fixed  four-layered  Goupillaud-type  strip  with  equal  layer  lengths  given  in 
(A.5),  can  be  extended  to  a  free-fixed  three-layered  non-Goupillaud-type  strip  with  wave  travel  time  ratios  1:2:1  by  choosing 
a2  =  1.  With  these  conditions,  (A.5)  becomes 


(1  +  oq)(l  +  a3)  cos2  -w  -  r3  cos  -w  +  (r3  -  (1  +  Gq)(l  +  a3)  +  2)  =  0,  (34) 

where  r3  =  (oqa3  —  1)  is  chosen  so  that  the  (cosine)  solutions  vary  in  the  interval  [—1,  1].  This  can  be  achieved  when 
r3  =  0,  which  implies  that  (1  +  Qq)(l  +  a3)  >  2.  As  a  result,  (34)  reduces  to: 


0  <  cos2 


(1  +  aq)(l  +  a3)  —  2 
(1  +  «i)(l  +a3) 


5.  Conclusions 

The  resonance  frequency  spectrum  was  derived  for  an  m-layered  Goupillaud-type  elastic  medium  that  was  subjected  to 
a  discrete,  time-harmonic  sinusoidal  forcing  function  at  one  end,  and  held  fixed  at  the  other  end.  Analytical  stress  solutions 
were  obtained  from  a  coupled  first-order  system  of  difference  equations  using  z-transform  methods;  key  elements  of  this 
approach  were  that  the  determinant  of  the  resulting  global  matrix  |Am  |  in  the  z-space  was  palindromic  polynomial  with  real 

I  m  I 

coefficients,  and  that  its  zeros  were  distinct  and  on  the  unit  circle.  We  showed  in  (19)  that  the  angles  90  =  tc  and 
corresponding  to  the  zeros  of  |Am  |  on  the  unit  circle  generate  the  resonance  frequency  spectrum  for  the  discrete  model.  The 
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0-angle  formulas  for  selected  designs  can  be  recovered  from  (10)  and  ( 12)— ( 14),  as  well  as  [20],  In  particular,  for  m  odd,  the 
positive  angles  coterminal  with  90  =  n  represent  resonance  frequencies  which  do  not  depend  on  the  material  properties. 
In  general,  the  0-angles,  hence  the  resonance  frequencies  for  the  discrete  model,  depend  only  on  (certain  combinations  of) 
the  impedance  ratios  and  not  on  any  other  parameters  such  as  the  length  of  the  strip,  et  cetera.  According  to  (19),  as  long  as 
m-layered  designs  with  different  impedance  ratios  have  their  0-angles  in  common,  they  are  guaranteed  to  have  the  same 
resonance  frequency  spectrum  for  the  discrete  model. 

The  natural  frequency  spectrum  of  a  free-fixed  m-layered  Goupillaud-type  strip  was  derived  from  the  resonance 
frequency  spectrum  of  the  discrete  model  (19)  using  the  linear  transformation  for  frequencies  (27).  An  alternative  derivation 
was  given  in  Section  3.1  by  analytically  solving  the  frequency  equation  with  up  to  five  layers.  Based  on  (28),  the  natural 
frequency  spectrum  of  a  free-fixed  m-layered  Goupillaud  type  strip  depends  on  (combinations  of)  the  impedance  ratios  and 
it  is  inversely  proportional  to  the  equal  wave  travel  time  for  each  layer  + .  Designs  with  these  quantities  in  common  share  the 
same  natural  frequency  spectrum  for  the  free-fixed  boundary  conditions,  see  Section  4. 1  for  illustrations.  The  transformation 
of  variables  (2)  was  successfully  applied  in  Sections  2  and  3,  to  simplify  the  problem  under  consideration  and  the  frequency 
equation. 

A  possible  design  modification  that  provides  a  desired  resonance  frequency  spectrum  for  the  discrete  model  was 
demonstrated  in  Section  4.2.  When  the  resonance  frequency  results  are  combined  with  the  stress  optimization  results 
obtained  in  [20],  we  conclude  that  for  a  given  optimal  design  of  a  Goupillaud-type  layered  strip  with  one  fixed  end, 
depending  on  the  type  of  loading  at  the  other  end,  one  can  get  either  the  best  or  the  worst  outcome  when  it  comes  to  the 
stress  amplitude.  Although  this  work  focused  on  a  Goupillaud-type  layered  elastic  strip,  the  frequency  results  were  shown  to 
extend  to  a  non-Goupillaud-type  layered  strip  with  integer  layer  length  ratios.  The  stress  formulas  involving  trigonometric 
sums  can  also  be  interpreted  using  the  Chebyshev  polynomials  of  the  first  and  second  kind. 

So  far  we  have  shown  that  the  palindromic  polynomial  with  real  coefficients  |Am|  has  all  its  zeros  on  the  unit  circle  for 
1  <  m  <  5,  see  Section  2.2.  Whether  this  is  true  for  any  in-layered  Goupillaud-type  strip  or  only  for  certain  classes  of 
designs  remains  to  be  investigated.  Preliminary  work  on  this  topic  involving  tridiagonal  Toeplitz  matrices  was  discussed  in 
Section  2.2.  In  addition,  the  experimental  validation  of  the  frequency  results  proposed  here  would  be  a  natural  and  valuable 
extension  of  this  work. 
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Appendix  A.  Solutions  of  the  frequency  equations  for  a  free-fixed  m-layered  Goupillaud-type  elastic  strip  (2  <  m  <  5) 

A.l.  Natural  frequency  spectrum  for  the  free-fixed  two-layer  case 

After  a  few  mathematical  manipulations  involving  trigonometric  identities,  the  frequency  equation  (33)  for  m  =  2 
becomes 

Qq  —  1  /i  —  2 

cosr&>= - =  - ,  where  xi  =  O  +  cq)  >  1  or  oq  >  0.  (A.l) 

ai  +  1  Xi 

This  suggests  the  following  natural  frequency  spectrum 


1 

-\  ( Xi  —  2\ 

COS  -  +  2l\7Z 

i 

,  a>i2  =  -  ■ 

-cos  1  (  — - —  J  +  2(/2  +  1)7T 

r 

L  V  xi  /  J 

T 

L  V  xi  /  J 

for  /i,  l2  =  0,  1,2, _ Here  oq  represents  the  impedance  ratio  between  the  first  and  the  second  layer,  while  xi  =  1  +  <*i- 

When  both  layers  are  occupied  by  the  same  material  (oq  =  1),  Eq.  (A.l)  recovers  a  known  literature  result  for  the  natural 
frequency  spectrum  of  a  homogeneous  strip  of  length  r  given  in  [10],  In  this  case,  the  frequency  equation  (A.l)  becomes 
cos  to)  =  0  which  implies  the  natural  frequency  spectrum  a>i  =  l2l+r>  ■  | ,  /  =  0,  1,  2  . . ..  Notice  that  (A.l)  implies  that 
cos  rw  7^  ±1. 


A.2.  Natural  frequency  spectrum  for  the  free-fixed  three-layer  case 


Using  the  fact  that  |Am|  =  |Am|,  the  frequency  equation  (33)  for  m  =  3  becomes  equivalent  to  solving  the  following 
equations: 


2r  2r  X2  —  2 

cos  — w  =  —  1  or  cos  — a)  = - 

3  3  X2 


where  x2  =  (1  +  oq)(l  +  oq)  >  1. 


(A.3) 
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Table  B.l 

Optimal  base  angle  values  for  the  m-layer  case. 


m 

2 

3 

4 

5 

Optimal  base  angles 

fll.opt  =  J 

n  _  tt 

V\,opt  —  2j — 1 

@i,opt 
@2,  opt 

_  7T 

-  2 r 
=*-§ 

$l,opf  — 
@2, opt  — 

71 

2/+1  ’ 

(2j— 1)tt 

2/+1 

From  (A.3)  the  frequency  spectrum  for  the  three-layer  is 


,  3tt 

=  (2 10  +  1)  — , 


a>q  = 


2r 


COS  1  (  — -  )  +  2/i7T 

X2 


COu  = 


2t 


-  cos  1  |  — — 2  )  +  2 (l2  +  1  )n 


where  l0,  l,,  l2  =  0,  1,  2, . . .. 

A.3.  Natural  frequency  spectrum  for  the  free-fixed  four-layer  case 

Use  of  the  double  angle  formula  for  the  cosine  for  the  four-layer  case  (m  =  4),  results  in 

2  T  2r3  T  2r3  —  /3  +  4 

cos  —co - cos  —a>  -I - =  0. 

2  Xs  2  xs 

For  any  given  values  of  the  design  parameters  /3  and  A.  the  frequency  spectrum  is: 


a>ij  = 


a>b  = 


cos 


/  A  ±  V(X3  -  A)2  -  4X3 


X3 


+ 


—  COS 


— i  /  r3  ±  V(X3  -  r3)2  -  4X3 


/i,/2  =  0,  1,2.... 


X3 


+  2(/2  +  1)jt 


(A.4) 


(A.5) 


(A.6) 


When  the  first  two  layers  are  occupied  by  one  material  (a 3  =  1)  and  the  other  two  are  occupied  by  another  material 
{oi2  1,  oi3  =  1),  the  four-layer  case  reduces  to  the  two-layer  case  and  (A.5)  reduces  to  (A.l).  From  (A.5)  if  follows  that 

2  r  a2 

cos  —co  = - , 

2  oc2  +  1 

which  after  applying  the  double  angle  formula  for  cosine  becomes  consistent  with  (A.l), 
of2  —  1 


COS  TO)  = 


a2  +  1 ' 


A.4.  Natural  frequency  spectrum  for  the  free-fixed  five-layer  case 

For  m  =  5,\A5\  =  cos  ™  ■  [X4  cos4  —  (X4  +  A)  c°s2  f  to  +  (A  +  1)]  =  0,  and  we  obtain 


2r 

COS  - ft)  : 

5 


2r 


2  A 


—  1  or  cos^  — co - 1  cos  — ft>  + 

5  X4  5 


2r  2  A  —  X4  "F  4 


0. 


X4 


(A.7) 


The  frequency  spectrum  can  then  be  obtained  for  any  given  values  of  X4  and  A.  similar  to  the  four  layer  case. 

Appendix  B.  Optimal  base  angles 

For  an  m-layered  Goupillaud-type  elastic  strip  subjected  to  a  Heaviside  loading  at  one  end  and  held  fixed  at  the  other 
end,  the  optimal  designs  have  the  smallest  stress  amplitude  of  double  the  loading.  These  optimal  designs  are  characterized 
by  the  optimal  values  in  Table  B.l  of  their  base  angles  for  j  =  2,3,4,...: 


Appendix  C.  Material  properties 

The  elastic  modulus  E  and  density  p  of  the  materials  used  in  Section  4.1  (see  [29])  are  as  follows: 

Aluminum  alloy:  E  =  70  GPa  and  p  =  2500  kg/m3;  Brass:  E  =  20  GPa  and  p  =  984.19  kg/m3;  Copper  alloy:  £  =  10  GPa 
and  p  =  515.19  kg/m3;  Lead:  £  =  14  GPa  and  p  =  11,  340  kg/m3;  Tungsten  alloy:  £  =  275  GPa  and  p  =  19,  610  kg/m3; 
Steel:  £  =  200  GPa  and  p  =  8000  kg/m3. 
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J  LAVERY 
POBOX  12211 

RESEARCH  TRIANGLE  PARK  NC 
27709-2211 

1  COMMANDER 

US  ARMY  RSRCH  OFC 
RDRL  ROE  M 
D  STEPP 
POBOX  12211 

RESEARCH  TRIANGLE  PARK  NC 
27709-221 

5  SOUTHWEST  RSRCH  INST 
C  ANDERSON 
K  DANNEMANN 
T  HOLMQUIST 
G  JOHNSON 
J  WALKER 
PO  DRAWER  28510 
SAN  ANTONIO  TX  78284 


1  ERDC 

US  ARMY  CORPS  OF  ENGINEERS 
USACEGSL 
P  PAP ADOS 
7701  TELEGRAPH  RD 
ALEXANDRIA  VA  22315 

1  AFOSR/NL 

875  NORTH  RANDOLPH  ST 
SUITE  325  RM  3112 
F  FAHROO 

ARLINGTON  VA  22203 

5  NAVAL  RESEARCH  LAB 
E  R  FRANCHI  CODE  7100 
M  H  ORR  CODE  7120 
J  A  BUCARO  CODE  7130 
J  S  PERKINS  CODE  7140 
S  A  CHIN  BING  CODE  7180 
4555  OVERLOOK  AVE  SW 
WASHINGTON  DC  20375 

1  UNIV  OF  MISSISSIPPI 
DEPT  OF  MECH  ENGRG 
A  M  RAJENDRAN 
201-B  CARRIER  HALL 
UNIVERSITY  MS  38677 

2  SRI  INTERNATIONAL 
D  CURRAN 

D  SHOCKEY 

333  RAVENS  WOOD  AVE 

MENLO  PARK  CA  94025 

1  VIRGINIA  POLYTECHNIC  INST 
COLLEGE  OF  ENGRG 
RBATRA 

BLACKSBURG  VA  24061-0219 

1  JOHNS  HOPKINS  UNIV 
DEPT  OF  MECH  ENGRG 
K  T  RAMESH 
LATROBE  122 
BALTIMORE  MD  21218 

1  INST  OF  ADVANCED  TECH 
UNIV  OF  TX  AUSTIN 
S  BLESS 

3925  W  BRAKER  LN  STE  400 
AUSTIN  TX  78759-5316 


NO.  OF 

COPIES  ORGANIZATION 

1  APPLIED  RSCH  ASSOCIATES 
D  E  GRADY 

4300  SAN  MATEO  BLVD  NE 
STE  A220 

ALBUQUERQUE  NM  871 10 

1  INTERNATIONAL  RSRCH 
ASSOC  INC 

D  L  ORPHAL  CAGE  06EXO 
4450  BLACK  AVE 
PLEASANTON  CA  94566 

2  WASHINGTON  ST  UNIV 
INST  OF  SHOCK  PHYSICS 
Y  M  GUPTA 

J  ASAY 

PULLMAN  WA  99164-2814 

1  UNIV  OF  DAYTON 

RSRCH  INST 
NSBRAR 

300  COLLEGE  PARK 
MS  SPC  1911 
DAYTON  OH  45469 

1  TEXAS  A&M  UNIV 
DEPT  OF  GEOLOGY  & 
GEOPHYSICS  MS  3115 
F  CHESTER 

COLLEGE  STATION  TX  778431 

1  UNIV  OF  SAN  DIEGO 

DEPT  OF  MATH  &  CMPTR  SCI 
A  VELO 

5998  ALCALA  PARK 
SAN  DIEGO  C A  92110 

1  NATIONAL  INST  OF 

STANDARDS  &  TECHLGY 
BLDG  &  FIRE  RSRCH  LAB 
J  MAIN 

100  BUREAU  DR  MS  861 1 
GAITHERSBURG  MD  20899-8611 

1  MIT 

DEPT  ARNTCS  ASTRNTCS 
R  RADOVITZKY 
77  MASSACHUSETTS  AVE 
CAMBRIDGE  MA  02139 

1  PENN  STATE  UNIV 

DEPT  OF  ENGRG  SCI  &  MECH 
F  COSTANZO 

UNIVERSITY  PARK  PA  168023 


NO.  OF 

COPIES  ORGANIZATION 

1  UNIV  OF  TEXAS-PAN  AMERICAN 
COLLEGE  OF  ENGRG 
&  COMPUTER  SCI 
D  H  ALLEN 

1201  WEST  UNIVERSITY  DR 
EDINBURG,  TX  78539-2999 

1  CLEMSON  UNIV 

DEPT  OF  MECH  ENGRG 
M  GRUJICIC 

241  ENGRG  INNOVATION  BLDG 
CLEMSON  SC  29634-0921 

1  UNIV  OF  DELAWARE 

DEPT  ELECTRICAL  &  CMPTR 
ENGRG 
D  WEILE 

NEWARK  DE  19716 

7  UNIV  OF  NEBRASKA 
DEPT  OF  ENGRG  MECH 
F  BOBARU 
Y  DZENIS 
G  GOGOS 
M  NEGAHBAN 
R  FENG 
J  TURNER 
Z  ZHANG 
LINCOLN  NE  68588 

1  WORCESTER  POLYTECHNIC  INST 
MATHEMATICAL  SCI 
K  LURIE 

WORCESTER  MA  01609 

4  UNIV  OF  UTAH 
DEPT  OF  MATH 
A  CHERKAEV 
E  CHERKAEV 
E  S  FOLIAS 
R  BRANNON 

SALT  LAKE  CITY  UT  841 12 

4  UNIV  OF  DELAWARE 
DEPT  OF  MECH  ENGRG 
T  BUCHANAN 
T  W  CHOU 
A  KARLSSON 
M  SANT ARE 
126  SPENCER  LAB 
NEWARK  DE  19716 


NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


1  NORTHWESTERN  UNIV 

DEPT  OF  CIVIL  &  ENVIRON  ENGRG 
Z  BAZANT 

2145  SHERIDAN  RD  A135 
EVANSTON  IL  60208-3109 

1  UNIV  OF  DELAWARE 

CTR  FOR  COMPST  MATRLS 
J  GILLESPIE 
NEWARK  DE  19716 

1  LOUISIANA  STATE  UNIV 
R  LIPTON 

304  LOCKETT  HALL 
BATON  ROUGE  LA  70803-4918 

1  UNIV  OF  ILLINOIS 

DEPT  OF  MECHL  SCI  &  ENGRG 
A  F  VAKAKIS 
1206  W  GREEN  ST  MC  244 
URBANA  CHAMPAIGN  IL  61801 

1  UNIV  OF  ILLINOIS 
ARSPC  ENGRG 
J  LAMBROS 

104  S  WRIGHT  ST  MC  236 
URBANA  CHAMPAIGN  IL  61801 

1  T  W  WRIGHT 

4906  WILMSLOW  RD 
BALTIMORE  MD  21210 

4  ADELPHI  LABORATORY  CTR 

C  CHABALOWSKI 
J  CHANG 
O  OCHOA 
R  SKAGGS 
BLDG  205 

2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 


ABERDEEN  PROVING  GROUND 

89  DIR  USARL 
RDRL  CIH  C 
J  CAZAMIAS 
P CHUNG 
D  GROVE 
J  KNAP 
RDRL  WM 
B  FORCH 
SKARNA 
J  MCCAULEY 


P  PLOSTINS 
RDRL  WML 
J  NEWILL 
M  ZOLTOSKI 
RDRL  WML  B 
I BATYREV 
S  IZVYEKOV 
B  RICE 

R  PESCE  RODRIGUEZ 
D  TAYLOR 
N  TRIVEDI 
N  WEINGARTEN 
RDRL  WML  D 
P  CONROY 
M  NUSCA 
RDRL  WML  E 
P  WEINACHT 
RDRL  WML  F 
D  LYON 
RDRL  WML  G 
M  BERMAN 
W  DRYSDALE 
RDRL  WML  H 
D  SCHEFFLER 
S  SCHRAML 
B  SCHUSTER 
RDRL  WMM 
J  BEATTY 
R  DOWDING 
J  ZABINSKI 
RDRL  WMM  A 
J  SANDS 
J  TZENG 
E  WETZEL 
RDRL  WMM  B 
T  BOGETT1 
B  CHEESEMAN 
C  FOUNTZOULAS 
G  GAZONAS 
D  HOPKINS 
R  KARKKAINEN 
PMOY 
B  POWERS 
C  RANDOW 
T  SANO 
F  TAVAZZA 
M  VANLANDINGHAM 
R  WILDMAN 
C  YEN 

RDRL  WMM  C 
J  LA  SCALA 
RDRL  WMM  D 
ECHIN 
KCHO 


NO.  OF 

COPIES  ORGANIZATION 


RDRL  WMM  E 
J  ADAMS 
M  COLE 
T  JESSEN 
J  LASALVIA 
P  PATEL 
J  SINGH 
RDRL  WMM  F 
L  KECSKES 
H  MAUPIN 
RDRL  WML  G 
J  ANDZELM 
A  RAWLETT 
RDRL  WMP 
P  BAKER 
S  SCHOENFELD 
RDRL  WMP  A 
B  RINGERS 
RDRL  WMP  B 
C  HOPPEL 
R  KRAFT 
S  SAT APATHY 
M  SCHEIDLER 
T  WEERASOORIY A 
RDRL  WMP  C 
R  BECKER 
S  BILYK 
T  BJERKE 
D  CASEM 
J  CLAYTON 
M  GREENFIELD 
B  LEAVY 
M  RAFTENBERG 
S  SEGLETES 
RDRL  WMP  D 
R  DONEY 
D  KLEPONIS 
J  RUNYEON 
B  SCOTT 
H  MEYER 
RDRL  WMP  E 
M  BURKINS 
B  LOVE 
RDRL  WMP  F 
A  FRYDMAN 
N  GNIAZDOWSKI 
R  GUPTA 
RDRL  WMP  G 
N  ELDREDGE 
D  KOOKER 
S  KUKUCK 


